Los datos provienen de la página web Kaggle. El dataset contiene información sobre setas comestibles y venenosas. El dataset contiene 8124 observaciones y 24 variables. Las variables con las que trabajaremos son las siguientes:
1. class: edible=e, poisonous=p
2. cap-diameter: continuous
3. cap-shape: bell=b, conical=c, convex=x, flat=f,
knobbed=k, sunken=s
4. cap-surface: fibrous=f, grooves=g, scaly=y,
smooth=s
5. cap-color: brown=n, buff=b, cinnamon=c, gray=g,
green=r, pink=p, purple=u, red=e, white=w, yellow=y
6. bruises: bruises=t, no=f
7. odor: almond=a, anise=l, creosote=c, fishy=y,
foul=f, musty=m, none=n, pungent=p, spicy=s
8. gill-attachment: attached=a, descending=d, free=f,
notched=n
9. gill-spacing: close=c, crowded=w, distant=d
10. gill-size: broad=b, narrow=n
11. gill-color: black=k, brown=n, buff=b, chocolate=h,
gray=g, green=r, orange=o, pink=p, purple=u, red=e, white=w,
yellow=y
12. stalk-shape: enlarging=e, tapering=t
13. stalk-root: bulbous=b, club=c, cup=u, equal=e,
rhizomorphs=z, rooted=r, missing=?
14. stalk-surface-above-ring: fibrous=f, scaly=y,
silky=k, smooth=s
15. stalk-surface-below-ring: fibrous=f, scaly=y,
silky=k, smooth=s
16. stalk-color-above-ring: brown=n, buff=b,
cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y
17. stalk-color-below-ring: brown=n, buff=b,
cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y
18. veil-type: partial=p, universal=u
19. veil-color: brown=n, orange=o, white=w,
yellow=y
20. ring-number: none=n, one=o, two=t
21. ring-type: cobwebby=c, evanescent=e, flaring=f,
large=l, none=n, pendant=p, sheathing=s, zone=z
22. spore-print-color: black=k, brown=n, buff=b,
chocolate=h, green=r, orange=o, purple=u, white=w, yellow=y
23. population: abundant=a, clustered=c, numerous=n,
scattered=s, several=v, solitary=y
24. habitat: grasses=g, leaves=l, meadows=m, paths=p,
urban=u, waste=w, woods=d
install.packages("caret")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/caret_6.0-93.tgz'
Content type 'application/x-gzip' length 3578325 bytes (3.4 MB)
==================================================
downloaded 3.4 MB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
install.packages("tidyverse")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/tidyverse_1.3.2.tgz'
Content type 'application/x-gzip' length 425892 bytes (415 KB)
==================================================
downloaded 415 KB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
install.packages("plotly")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/plotly_4.10.1.tgz'
Content type 'application/x-gzip' length 3179590 bytes (3.0 MB)
==================================================
downloaded 3.0 MB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
install.packages("dplyr")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/dplyr_1.0.10.tgz'
Content type 'application/x-gzip' length 1327807 bytes (1.3 MB)
==================================================
downloaded 1.3 MB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
install.packages("factoextra")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/factoextra_1.0.7.tgz'
Content type 'application/x-gzip' length 415155 bytes (405 KB)
==================================================
downloaded 405 KB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
install.packages("dendextend")
trying URL 'https://cran.rstudio.com/bin/macosx/big-sur-arm64/contrib/4.2/dendextend_1.16.0.tgz'
Content type 'application/x-gzip' length 3895351 bytes (3.7 MB)
==================================================
downloaded 3.7 MB
The downloaded binary packages are in
/var/folders/v0/47swnvh93cl0_bw1dw8tl_0h0000gn/T//RtmpBAcuFX/downloaded_packages
library(caret)
Loading required package: ggplot2
Loading required package: lattice
Registered S3 method overwritten by 'data.table':
method from
print.data.table
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ──────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2
✔ purrr 0.3.5 ── Conflicts ─────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ purrr::lift() masks caret::lift()
library(plotly)
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: ‘plotly’
The following object is masked from ‘package:ggplot2’:
last_plot
The following object is masked from ‘package:stats’:
filter
The following object is masked from ‘package:graphics’:
layout
library(dplyr)
library(cluster)
library(factoextra)
Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(dendextend)
---------------------
Welcome to dendextend version 1.16.0
Type citation('dendextend') for how to cite the package.
Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/
Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
You may ask questions at stackoverflow, use the r and dendextend tags:
https://stackoverflow.com/questions/tagged/dendextend
To suppress this message use: suppressPackageStartupMessages(library(dendextend))
---------------------
Attaching package: ‘dendextend’
The following object is masked from ‘package:stats’:
cutree
En primer lugar, realizamos la lectura del fichero csv, separándolo por “;” y mostramos las primeras 6 filas del fichero.
mushroom <- read.csv("./data/data.csv", sep = ";")
head(mushroom)
Echamos un vistazo a las características de los atributos del dataset. En el caso de las variables numéricas, se puede observar valores como el mínimo, máximo, media, desviación estándar, etc.
summary(mushroom)
class cap.diameter cap.shape cap.surface
Length:61069 Min. : 0.380 Length:61069 Length:61069
Class :character 1st Qu.: 3.480 Class :character Class :character
Mode :character Median : 5.860 Mode :character Mode :character
Mean : 6.734
3rd Qu.: 8.540
Max. :62.340
cap.color does.bruise.or.bleed gill.attachment gill.spacing
Length:61069 Length:61069 Length:61069 Length:61069
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
gill.color stem.height stem.width stem.root stem.surface
Length:61069 Min. : 0.000 Min. : 0.00 Length:61069 Length:61069
Class :character 1st Qu.: 4.640 1st Qu.: 5.21 Class :character Class :character
Mode :character Median : 5.950 Median : 10.19 Mode :character Mode :character
Mean : 6.582 Mean : 12.15
3rd Qu.: 7.740 3rd Qu.: 16.57
Max. :33.920 Max. :103.91
stem.color veil.type veil.color has.ring
Length:61069 Length:61069 Length:61069 Length:61069
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
ring.type spore.print.color habitat season
Length:61069 Length:61069 Length:61069 Length:61069
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
Debido a que el dataset contiene valores vacíos, los sustituiremos por valores NA para poder trabajar con ellos.
mushroom[mushroom == ""] <- NA
Comprobamos que se han sustituido correctamente los valores vacíos por NA.
str(mushroom)
'data.frame': 61069 obs. of 21 variables:
$ class : chr "p" "p" "p" "p" ...
$ cap.diameter : num 15.3 16.6 14.1 14.2 14.6 ...
$ cap.shape : chr "x" "x" "x" "f" ...
$ cap.surface : chr "g" "g" "g" "h" ...
$ cap.color : chr "o" "o" "o" "e" ...
$ does.bruise.or.bleed: chr "f" "f" "f" "f" ...
$ gill.attachment : chr "e" "e" "e" "e" ...
$ gill.spacing : chr NA NA NA NA ...
$ gill.color : chr "w" "w" "w" "w" ...
$ stem.height : num 16.9 18 17.8 15.8 16.5 ...
$ stem.width : num 17.1 18.2 17.7 16 17.2 ...
$ stem.root : chr "s" "s" "s" "s" ...
$ stem.surface : chr "y" "y" "y" "y" ...
$ stem.color : chr "w" "w" "w" "w" ...
$ veil.type : chr "u" "u" "u" "u" ...
$ veil.color : chr "w" "w" "w" "w" ...
$ has.ring : chr "t" "t" "t" "t" ...
$ ring.type : chr "g" "g" "g" "p" ...
$ spore.print.color : chr NA NA NA NA ...
$ habitat : chr "d" "d" "d" "d" ...
$ season : chr "w" "u" "w" "w" ...
Comprobaremos la cantidad de valores nulos que hay en cada variable. Para ello, utilizaremos la función colSums(is.na(mushroom)), que nos devolverá la suma de valores nulos de cada variable. Esta información nos ayudará a decidir si será conveniente eliminar dichas variables o no. Podemos observar que existen 5 variables donde más del 50% de los valores son nulos. Estas son: stem.surface, veil.color, spore.print.color, stem.root, veil.type. Esta información nos ayudará a decidir si será conveniente eliminar dichas variables o no. En este caso, dichas variables serán eliminadas ya que, en caso de imputar los valores nulos, se perdería mucha información.
colSums(is.na(mushroom))
class cap.diameter cap.shape cap.surface
0 0 0 14120
cap.color does.bruise.or.bleed gill.attachment gill.spacing
0 0 9884 25063
gill.color stem.height stem.width stem.root
0 0 0 51538
stem.surface stem.color veil.type veil.color
38124 0 57892 53656
has.ring ring.type spore.print.color habitat
0 2471 54715 0
season
0
En este caso, decidimos eliminar aquellas variables con más del 50% de valores nulos ya que, en caso de imputar los valores nulos, se perdería mucha información. Para ello, en primer lugar, obtendremos el nombre de las columnas que queremos eliminar.
nacols <- colnames(mushroom)[colSums(is.na(mushroom)) > nrow(mushroom) / 2]
print(nacols)
[1] "stem.root" "stem.surface" "veil.type" "veil.color"
[5] "spore.print.color"
Eliminamos las variables no deseadas.
mushroom <- mushroom[, !names(mushroom) %in% nacols]
Comprobamos que se han eliminado las variables con más del 50% de valores nulos.
print(colnames(mushroom))
[1] "class" "cap.diameter" "cap.shape"
[4] "cap.surface" "cap.color" "does.bruise.or.bleed"
[7] "gill.attachment" "gill.spacing" "gill.color"
[10] "stem.height" "stem.width" "stem.color"
[13] "has.ring" "ring.type" "habitat"
[16] "season"
A continuación, separamos las variables numéricas de las categóricas.
colsnames <- colnames(mushroom)
numerical_features <- c("cap.diameter", "stem.height", "stem.width")
categorical_features <- colsnames[!colsnames %in% numerical_features]
print(categorical_features)
[1] "class" "cap.shape" "cap.surface"
[4] "cap.color" "does.bruise.or.bleed" "gill.attachment"
[7] "gill.spacing" "gill.color" "stem.color"
[10] "has.ring" "ring.type" "habitat"
[13] "season"
print(numerical_features)
[1] "cap.diameter" "stem.height" "stem.width"
Visualizamos la distribución de las variables categóricas a través de histogramas. Observaremos los posibles valores de cada variable categórica junto con su frecuencia de aparición en el dataset.
for (i in categorical_features) {
print(ggplot(mushroom, aes_string(x = i)) +
geom_bar(fill = "#7fd6d9") +
geom_text(stat = "count", aes(label = scales::percent(..count.. / nrow(mushroom)), vjust = -0.25)) +
labs(x = i, y = "Percentage") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)))
}
Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
Please use tidy evaluation ideoms with `aes()`
Tras visualizar los histogramas de cada variable, llegamos a las siguientes conclusiones:
A continuación, visualizaremos la distribución de las variables numéricas a través de una gráfica generada con “featurePlot”. Esta requiere que la variable objetivo sea de tipo factor, por lo que hacemos la conversión.
mushroom$class <- as.factor(mushroom$class)
featurePlot(x = mushroom[, numerical_features], y = mushroom$class, plot = "strip")
Con la gráfica anterior, se puede observar que para valores altos de cap.diameter, stem.height y stem.width, la probabilidad de que la clase sea “e” es mayor que la de “p”.
Para decidir si es conveniente eliminar alguna variable, se puede utilizar la función “nearZeroVar”. Esta función devuelve un vector con los índices de las variables que tienen una varianza cercana a 0.
near_zero_col <- nearZeroVar(mushroom, saveMetrics = FALSE)
colnames(mushroom)[c(near_zero_col)]
[1] "ring.type"
Como conclusión, la variable “ring.type” podría ser eliminada ya que tiene una varianza cercana a 0. A continuación, visualizaremos la correlación existente con la variable dependiente “class”.
print(ggplot(mushroom, aes_string(x = "ring.type")) +
geom_bar(aes(fill = class)))
Tras visualizar la gráfica anterior, se puede observar que la distribución de la variable dependiente “class” es similar en casi todos los valores de “ring.type”. Sin embargo, en el caso de “ring.type” = “z” y “ring.type” = “m”, la distribución de la variable dependiente “class” es diferente. Por lo tanto, se decide no eliminar la variable “ring.type”.
Como se ha comentado anteriormente, los valores nulos de las variables categóricas se imputarán a través de la moda de cada variable. Esto es debido a que la función preProcess() de la librería caret no permite imputar valores nulos de variables categóricas.
for (i in categorical_features) {
mushroom[, i][is.na(mushroom[, i])] <- names(which.max(table(mushroom[, i])))
}
Podemos comprobar que ya no existen valores nulos en las variables categóricas.
colSums(is.na(mushroom))
class cap.diameter cap.shape cap.surface
0 0 0 0
cap.color does.bruise.or.bleed gill.attachment gill.spacing
0 0 0 0
gill.color stem.height stem.width stem.color
0 0 0 0
has.ring ring.type habitat season
0 0 0 0
Para escalar las variables numéricas, se utilizará la función preProcess() de la librería caret. Esta función devuelve un objeto de tipo “preProcess” que contiene la información necesaria para escalar las variables numéricas. A continuación, se escalarán las variables numéricas y se eliminarán las variables originales.
range_numeric <- preProcess(mushroom[, numerical_features], method = c("range"))
mushroom[, numerical_features] <- predict(range_numeric, newdata = mushroom[, numerical_features])
str(mushroom)
'data.frame': 61069 obs. of 16 variables:
$ class : Factor w/ 2 levels "e","p": 2 2 2 2 2 2 2 2 2 2 ...
$ cap.diameter : num 0.24 0.262 0.221 0.223 0.23 ...
$ cap.shape : chr "x" "x" "x" "f" ...
$ cap.surface : chr "g" "g" "g" "h" ...
$ cap.color : chr "o" "o" "o" "e" ...
$ does.bruise.or.bleed: chr "f" "f" "f" "f" ...
$ gill.attachment : chr "e" "e" "e" "e" ...
$ gill.spacing : chr "c" "c" "c" "c" ...
$ gill.color : chr "w" "w" "w" "w" ...
$ stem.height : num 0.5 0.53 0.525 0.465 0.487 ...
$ stem.width : num 0.164 0.175 0.171 0.154 0.166 ...
$ stem.color : chr "w" "w" "w" "w" ...
$ has.ring : chr "t" "t" "t" "t" ...
$ ring.type : chr "g" "g" "g" "p" ...
$ habitat : chr "d" "d" "d" "d" ...
$ season : chr "w" "u" "w" "w" ...
Primero se dividirá el dataset en dos conjuntos: uno de entrenamiento y otro de test. El conjunto de entrenamiento se utilizará para entrenar los modelos y el conjunto de test se utilizará para evaluarlos.
library(caTools)
set.seed(18)
split <- sample.split(mushroom$class, SplitRatio = 0.8)
training_set <- subset(mushroom, split == TRUE)
test_set <- subset(mushroom, split == FALSE)
table(training_set$class)
e p
21745 27110
table(test_set$class)
e p
5436 6778
(Redactar) La regresión logística permite predecir el resultado de una variable categórica en función de las variables independientes o predictoras.
rl_classiffier <- glm(class ~ ., family = binomial, data = training_set)
summary(rl_classiffier)
Call:
glm(formula = class ~ ., family = binomial, data = training_set)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.98113 -0.77757 0.00045 0.75561 2.99676
Coefficients: (2 not defined because of singularities)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -16.146622 310.494270 -0.052 0.958526
cap.diameter -3.097299 0.273141 -11.340 < 2e-16 ***
cap.shapec -1.587320 0.089853 -17.666 < 2e-16 ***
cap.shapef -1.548574 0.056840 -27.244 < 2e-16 ***
cap.shapeo -0.505473 0.102060 -4.953 7.32e-07 ***
cap.shapep -1.355581 0.077525 -17.486 < 2e-16 ***
cap.shapes -1.914206 0.067428 -28.389 < 2e-16 ***
cap.shapex -1.587070 0.052730 -30.098 < 2e-16 ***
cap.surfacee 0.777609 0.081079 9.591 < 2e-16 ***
cap.surfaceg -0.469069 0.067815 -6.917 4.62e-12 ***
cap.surfaceh -0.869361 0.063758 -13.635 < 2e-16 ***
cap.surfacei 1.878605 0.116592 16.113 < 2e-16 ***
cap.surfacek 3.205155 0.106739 30.028 < 2e-16 ***
cap.surfacel -1.338110 0.100947 -13.256 < 2e-16 ***
cap.surfaces -0.913831 0.059362 -15.394 < 2e-16 ***
cap.surfacet -0.037125 0.050461 -0.736 0.461905
cap.surfacew -0.483820 0.084640 -5.716 1.09e-08 ***
cap.surfacey -0.579375 0.064113 -9.037 < 2e-16 ***
cap.colore 1.986990 0.109951 18.072 < 2e-16 ***
cap.colorg 0.582318 0.108490 5.367 7.98e-08 ***
cap.colork 1.553011 0.131690 11.793 < 2e-16 ***
cap.colorl -0.129820 0.140094 -0.927 0.354103
cap.colorn 0.374666 0.099619 3.761 0.000169 ***
cap.coloro 1.627975 0.112745 14.439 < 2e-16 ***
cap.colorp 1.051060 0.127472 8.245 < 2e-16 ***
cap.colorr 2.925006 0.136340 21.454 < 2e-16 ***
cap.coloru 1.538956 0.121348 12.682 < 2e-16 ***
cap.colorw 0.861936 0.104891 8.217 < 2e-16 ***
cap.colory 0.615058 0.102936 5.975 2.30e-09 ***
does.bruise.or.bleedt -0.147991 0.038250 -3.869 0.000109 ***
gill.attachmentd 0.644857 0.046979 13.726 < 2e-16 ***
gill.attachmente -0.922859 0.055651 -16.583 < 2e-16 ***
gill.attachmentf 0.816396 0.136750 5.970 2.37e-09 ***
gill.attachmentp -2.450165 0.060717 -40.354 < 2e-16 ***
gill.attachments 0.142043 0.052091 2.727 0.006395 **
gill.attachmentx 0.089014 0.043392 2.051 0.040231 *
gill.spacingd -0.425196 0.036863 -11.534 < 2e-16 ***
gill.spacingf NA NA NA NA
gill.colore 2.304077 0.156835 14.691 < 2e-16 ***
gill.colorf NA NA NA NA
gill.colorg 0.833125 0.124899 6.670 2.55e-11 ***
gill.colork 0.854961 0.136405 6.268 3.66e-10 ***
gill.colorn 1.350254 0.119303 11.318 < 2e-16 ***
gill.coloro 1.131719 0.123598 9.156 < 2e-16 ***
gill.colorp 0.848246 0.122841 6.905 5.01e-12 ***
gill.colorr 0.896493 0.144108 6.221 4.94e-10 ***
gill.coloru 1.319431 0.145603 9.062 < 2e-16 ***
gill.colorw 0.819802 0.115149 7.120 1.08e-12 ***
gill.colory 1.818426 0.119373 15.233 < 2e-16 ***
stem.height 3.580359 0.196103 18.258 < 2e-16 ***
stem.width -0.455339 0.225084 -2.023 0.043076 *
stem.colore 17.897727 310.494233 0.058 0.954033
stem.colorf 35.399323 334.067154 0.106 0.915610
stem.colorg 15.906809 310.494228 0.051 0.959142
stem.colork 19.545176 310.494262 0.063 0.949807
stem.colorl 15.591779 310.494286 0.050 0.959950
stem.colorn 17.491009 310.494224 0.056 0.955077
stem.coloro 16.096825 310.494231 0.052 0.958654
stem.colorp 19.131727 310.494251 0.062 0.950868
stem.colorr 17.537012 310.494253 0.056 0.954959
stem.coloru 17.419638 310.494223 0.056 0.955260
stem.colorw 16.082601 310.494224 0.052 0.958691
stem.colory 17.437089 310.494226 0.056 0.955215
has.ringt -0.006793 0.050561 -0.134 0.893124
ring.typef -0.911592 0.084553 -10.781 < 2e-16 ***
ring.typeg -0.499679 0.101534 -4.921 8.60e-07 ***
ring.typel -0.167323 0.100958 -1.657 0.097449 .
ring.typem -20.187930 230.159347 -0.088 0.930105
ring.typep 0.894990 0.103503 8.647 < 2e-16 ***
ring.typer -0.492646 0.105898 -4.652 3.29e-06 ***
ring.typez 16.741684 78.753199 0.213 0.831651
habitatg 0.548824 0.041040 13.373 < 2e-16 ***
habitath 0.199731 0.065306 3.058 0.002225 **
habitatl -0.501569 0.054900 -9.136 < 2e-16 ***
habitatm 0.103008 0.067655 1.523 0.127870
habitatp 17.002090 230.326776 0.074 0.941156
habitatu -17.054326 378.536784 -0.045 0.964065
habitatw -17.436003 215.700665 -0.081 0.935574
seasons -1.618439 0.069225 -23.379 < 2e-16 ***
seasonu 0.181688 0.025335 7.171 7.42e-13 ***
seasonw -1.277528 0.049110 -26.014 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 67137 on 48854 degrees of freedom
Residual deviance: 45160 on 48776 degrees of freedom
AIC: 45318
Number of Fisher Scoring iterations: 16
Una vez hemos aplicado la función glm(), obtenemos los valores residuales del modelo y los coeficientes de ajuste para cada una de las variables independientes. Además, obtenemos el p-value correspondiente a cada una de ellas. Las variables con dos o tres asteriscos aportan bastante relevancia como predictores; sin embargo, las variables que no poseen ninguno o uno son de menor relevancia.
A continuación, procedemos a predecir las clases del conjunto de entrenamiento y de validación. Tomamos como umbral 0,5, de forma que si la probabilidad queda por encima de dicho umbral es que es comestible, pero si queda por debajo es que es venenoso. Por último crearemos la matriz de confusión para valorar los resultados de las predicciones.
pred_train <- predict(rl_classiffier, newdata = training_set, type = "response")
Warning: prediction from a rank-deficient fit may be misleading
pred_train <- ifelse(pred_train > 0.5, "p", "e")
pred_train <- factor(pred_train, levels = c("e", "p"), labels = c("e", "p"))
confusion_m <- table(training_set$class, pred_train)
confusion_m
pred_train
e p
e 16489 5256
p 5573 21537
accuracy <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy
[1] 0.7783441
Una vez hemos visto el resultado que muestra la matriz de confusión, se observa que la predicción tiene una precisión del 77,8%, lo que puede resultar algo baja.
Repetimos el mismo proceso para el conjunto de test.
pred_test <- predict(rl_classiffier, newdata = test_set, type = "response")
Warning: prediction from a rank-deficient fit may be misleading
pred_test <- ifelse(pred_test > 0.5, "p", "e")
pred_test <- factor(pred_test, levels = c("e", "p"), labels = c("e", "p"))
confusion_m <- table(test_set$class, pred_test)
confusion_m
pred_test
e p
e 4087 1349
p 1420 5358
accuracy_rl <- sum(diag(confusion_m)) / sum(confusion_m)
accuracy_rl
[1] 0.7732929
Obtenemos una precisión con gran similitud a la del conjunto de training, un 77,3%. Para finalizar la regresión logística, vamos a graficar la curva ROC, la cual nos aporta la visualización de la relación entre los falsos y verdaderos positivos.
library(ROCR)
pred_rl_roc <- prediction(as.numeric(pred_test), as.numeric(test_set$class))
perf_rl_roc <- performance(pred_rl_roc, "tpr", "fpr")
perf_rl_auc <- performance(pred_rl_roc, "auc")
print(perf_rl_auc@y.values[[1]])
[1] 0.7711691
plot(perf_rl_roc)
Observando la curva ROC resultante podemos comentar que se mantiene por encima de la diagonal, lo que es buena señal, pero se aproxima a ella, pudiendo haber proporcionado resultados mejores.
Antes de nada, para poder aplicar knn, debemos transformar las variables categóricas en numéricas.
mushroom_num <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom_num <- as.data.frame(mushroom_num)
Visualizamos las variables categóricas codificadas y las dimensiones del dataset.
dim_mushroom <- dim(mushroom_num)
print(dim_mushroom)
[1] 61069 81
str(mushroom_num)
'data.frame': 61069 obs. of 81 variables:
$ class.p : num 1 1 1 1 1 1 1 1 1 1 ...
$ cap.diameter : num 0.24 0.262 0.221 0.223 0.23 ...
$ cap.shapec : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.shapef : num 0 0 0 1 0 0 1 0 1 1 ...
$ cap.shapeo : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.shapep : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.shapes : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.shapex : num 1 1 1 0 1 1 0 1 0 0 ...
$ cap.surfacee : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfaceg : num 1 1 1 0 0 1 0 0 1 1 ...
$ cap.surfaceh : num 0 0 0 1 1 0 1 1 0 0 ...
$ cap.surfacei : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfacek : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfacel : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfaces : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfacet : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfacew : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.surfacey : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colore : num 0 0 0 1 0 0 0 1 0 1 ...
$ cap.colorg : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colork : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colorl : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colorn : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.coloro : num 1 1 1 0 1 1 1 0 1 0 ...
$ cap.colorp : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colorr : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.coloru : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colorw : num 0 0 0 0 0 0 0 0 0 0 ...
$ cap.colory : num 0 0 0 0 0 0 0 0 0 0 ...
$ does.bruise.or.bleedt: num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.attachmentd : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.attachmente : num 1 1 1 1 1 1 1 1 1 1 ...
$ gill.attachmentf : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.attachmentp : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.attachments : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.attachmentx : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.spacingd : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.spacingf : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colore : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorf : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorg : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colork : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorn : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.coloro : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorp : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorr : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.coloru : num 0 0 0 0 0 0 0 0 0 0 ...
$ gill.colorw : num 1 1 1 1 1 1 1 1 1 1 ...
$ gill.colory : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.height : num 0.5 0.53 0.525 0.465 0.487 ...
$ stem.width : num 0.164 0.175 0.171 0.154 0.166 ...
$ stem.colore : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorf : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorg : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colork : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorl : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorn : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.coloro : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorp : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorr : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.coloru : num 0 0 0 0 0 0 0 0 0 0 ...
$ stem.colorw : num 1 1 1 1 1 1 1 1 1 1 ...
$ stem.colory : num 0 0 0 0 0 0 0 0 0 0 ...
$ has.ringt : num 1 1 1 1 1 1 1 1 1 1 ...
$ ring.typef : num 0 0 0 0 0 0 0 0 0 0 ...
$ ring.typeg : num 1 1 1 0 0 0 1 0 0 0 ...
$ ring.typel : num 0 0 0 0 0 0 0 0 0 0 ...
$ ring.typem : num 0 0 0 0 0 0 0 0 0 0 ...
$ ring.typep : num 0 0 0 1 1 1 0 1 1 1 ...
$ ring.typer : num 0 0 0 0 0 0 0 0 0 0 ...
$ ring.typez : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatg : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitath : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatl : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatm : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatp : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatu : num 0 0 0 0 0 0 0 0 0 0 ...
$ habitatw : num 0 0 0 0 0 0 0 0 0 0 ...
$ seasons : num 0 0 0 0 0 0 0 0 0 0 ...
$ seasonu : num 0 1 0 0 0 1 0 1 0 0 ...
$ seasonw : num 1 0 1 1 1 0 1 0 0 1 ...
library(caTools)
set.seed(18)
split <- sample.split(mushroom_num$class, SplitRatio = 0.8)
training_set_num <- subset(mushroom_num, split == TRUE)
test_set_num <- subset(mushroom_num, split == FALSE)
table(training_set_num$class)
0 1
21745 27110
table(test_set_num$class)
0 1
5436 6778
Para comenzar a aplicar el clasificador del vecino más cercano, haremos uso de la función knn() tomando como valor de k=5.
library(class)
set.seed(18)
pred_knn <- knn(train = training_set_num[, -1], test = test_set_num[, -1], cl = training_set_num$class, k = 5)
Una vez tenemos los resultados de las predicciones, cosntruimos la matriz de confusión.
summary(pred_knn)
0 1
5429 6785
confusionM <- table(test_set_num$class, pred_knn)
confusionM
pred_knn
0 1
0 5429 7
1 0 6778
accuracy_knn <- sum(diag(confusionM)) / sum(confusionM)
accuracy_knn
[1] 0.9994269
Como resultado de la matriz de confusión, tenemos una precisión del 99%, resultando ser mejor que en la regresión logística.
library(ROCR)
pred_knn_roc <- prediction(as.numeric(pred_knn), as.numeric(test_set_num$class))
perf_knn_roc <- performance(pred_knn_roc, "tpr", "fpr")
perf_knn_auc <- performance(pred_knn_roc, "auc")
print(perf_knn_auc@y.values[[1]])
[1] 0.9993561
plot(perf_knn_roc)
Una vez obtenemos la curva ROC y el resultado del área de debajo de la misma, observamos que el resultado es muy positivo, ya que gráficamente se aproxima a la vertical.
Comenzamos aplicando el árbol de decisión mediante la función rpart() como se muestra a continuación.
library(rpart)
Attaching package: ‘rpart’
The following object is masked from ‘package:dendextend’:
prune
set.seed(18)
dt_classiffier <- rpart(class ~ ., data = training_set)
library(rattle)
Loading required package: bitops
Rattle: A free graphical interface for data science with R.
Version 5.5.1 Copyright (c) 2006-2021 Togaware Pty Ltd.
Type 'rattle()' to shake, rattle, and roll your data.
tree2 <- rpart(class ~ ., training_set,
method = "class",
control = rpart.control(cp = 0.00001)
)
fancyRpartPlot(tree2)
pruned <- prune(tree2, cp = 0.01)
fancyRpartPlot(pruned)
Construimos la matriz de confusión al igual que en los casos anteriores.
pred_dt <- predict(dt_classiffier, newdata = test_set, type = "class")
confusionM <- table(test_set$class, pred_dt)
confusionM
pred_dt
e p
e 4257 1179
p 804 5974
accuracy_dt <- sum(diag(confusionM)) / sum(confusionM)
accuracy_dt
[1] 0.8376453
Una vez tenemos el resultado para el valor de precisión en la predicción una vez aplicado el árbol de decisión, 83%, pasamos a aplicar de forma visual la construcción de la gráfica de la curva ROC y el cálculo del AUC. En este caso, continúa siendo mejor resultado que el obtenido mediante la regresión logística, pero quedando por debajo que el resultado obtenido al aplicar KNN.
library(ROCR)
pred_dt_roc <- prediction(as.numeric(pred_dt), as.numeric(test_set$class))
perf_dt_roc <- performance(pred_dt_roc, "tpr", "fpr")
perf_dt_auc <- performance(pred_dt_roc, "auc")
print(perf_dt_auc@y.values[[1]])
[1] 0.8322468
plot(perf_dt_roc)
Comenzamos aplicando el clasificador de Random Forest mediante la función llamada de la misma forma; es decir, randomForest(). En esta función, el valor correspondiente al parámetro ntree indica la cantidad de árboles de decisión que formarán parte del clasificador.
library(randomForest)
randomForest 4.7-1.1
Type rfNews() to see new features/changes/bug fixes.
Attaching package: ‘randomForest’
The following object is masked from ‘package:rattle’:
importance
The following object is masked from ‘package:dplyr’:
combine
The following object is masked from ‘package:ggplot2’:
margin
set.seed(18)
rf_classiffier <- randomForest(class ~ ., data = training_set, ntree = 250)
Realizamos las predicciones sobre el conjunto de datos y construimos la matriz de confusión.
pred_rf <- predict(rf_classiffier, newdata = test_set, type = "class")
confusionM <- table(test_set$class, pred_rf)
confusionM
pred_rf
e p
e 5431 5
p 0 6778
accuracy_rf <- sum(diag(confusionM)) / sum(confusionM)
accuracy_rf
[1] 0.9995906
De forma similar al clasificador KNN, obtenemos un valor de precisión del 99%. Definimos la curva ROC para este caso y realizamos el cálculo del área bajo la curva, los cuales nos simbolizan su altísima precisión.
library(ROCR)
pred_rf_roc <- prediction(as.numeric(pred_rf), as.numeric(test_set$class))
perf_rf_roc <- performance(pred_rf_roc, "tpr", "fpr")
perf_rf_auc <- performance(pred_rf_roc, "auc")
print(perf_rf_auc@y.values[[1]])
[1] 0.9995401
plot(perf_rf_roc)
Comenzamos a aplicar el clasificador Máquina de Soporte Vectorial haciendo uso de la función svm(). En esta función, los valores de los parámetros type y kernel hacen referencia al tipo de clasificador; es decir, que el kernel es de tipo radial y gaussiano.
library(e1071)
set.seed(18)
svm_classiffier <- svm(class ~ .,
data = training_set,
type = "C-classification", kernel = "radial"
)
La predicción y la matríz de confusión son entonces las mostradas a continuación.
pred_svm <- predict(svm_classiffier, newdata = test_set, type = "class")
confusionM <- table(test_set$class, pred_svm)
confusionM
pred_svm
e p
e 5129 307
p 240 6538
accuracy_svm <- sum(diag(confusionM)) / sum(confusionM)
accuracy_svm
[1] 0.9552153
Como se puede observar el valor de la predicción en este caso resulta ser del 95%. Para finalizar, construimos la curva ROC correspondiente en este caso
library(ROCR)
pred_svm_roc <- prediction(as.numeric(pred_svm), as.numeric(test_set$class))
perf_svm_roc <- performance(pred_svm_roc, "tpr", "fpr")
perf_svm_auc <- performance(pred_svm_roc, "auc")
print(perf_svm_auc@y.values[[1]])
[1] 0.954058
plot(perf_svm_roc)
Una vez aplicados los cinco distintos métodos de clasificación sobre nuestro dataset tras realizar el preprocesamiento, podemos concluir diciendo que los clasificadores de KNN y Random Forest son los que poseen mayor precisión y, por tanto, menor error de predicción. Sin embargo, con los resultados obtenidos por parte de ambos podemos detectar un posible sobreajuste, no beneficiando al modelo. Por esta razón, creemos que resulta más beneficioso sacrificar parte del valor de precisión, teniendo en cuenta algunos valores de falsos positivos y falsos negativos, como ocurre en el caso del clasificador de la Máquina de Soporte Vectorial, aprovechando así su capacidad de mayor generalización.
accuracy_comp <- matrix(c(accuracy_rl, accuracy_knn, accuracy_dt, accuracy_rf, accuracy_svm), ncol = 5)
barplot(accuracy_comp)
barplot(accuracy_comp,
main = "Accuracy Comparison",
xlab = "Accuracy (%)",
ylab = "Method",
names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
col = "darkred"
)
perf_auc <- matrix(c(perf_rl_auc@y.values[[1]], perf_knn_auc@y.values[[1]], perf_dt_auc@y.values[[1]], perf_rf_auc@y.values[[1]], perf_svm_auc@y.values[[1]]), ncol = 5)
barplot(perf_auc)
barplot(perf_auc,
main = "AUC Comparison",
xlab = "AUC (%)",
ylab = "Method",
names.arg = c("RL", "K-NN", "DT", "RF", "SVM"),
col = "darkred"
)
En este apartado se analizará el dataset a través de algoritmos de aprendizaje no supervisado. En concreto, se probarán los algoritmos k-means y clustering jerárquico. Para ambos algoritmos, se seguirá el siguiente esquema:
Para poder trabajar con algoritmos no supervisados será necesario que las variables sean numéricas. Para ello, se eliminarán las variables categóricas del dataset.
numerical_columns <- mushroom[, numerical_features]
En primer lugar, representaremos de forma gráfica la distribución inicial de los datos.
mushroom <- dummyVars(" ~ .", data = mushroom, fullRank = TRUE) %>% predict(mushroom)
mushroom <- as.data.frame(mushroom)
df <- as.data.frame(numerical_columns)
plot_ly(df,
x = ~cap.diameter, y = ~stem.height,
z = ~stem.width
) %>%
add_markers(size = 1.5)
Para aplicar el algoritmo kmeans se utilizará la función kmeans() de la librería cluster. Será necesario determinar el número de clusters óptimo. Para ello, se utilizará la función “kmeans” con diferentes valores de “centers” y se calculará la suma de cuadrados internos (within groups sum of squares) para cada valor de “centers”. A continuación, se representará la suma de cuadrados internos vs. número de clusters.
wss_per_k <- 0
for (i in 1:10) {
kmeans_aux <- kmeans(numerical_columns, center = i, nstar = 20)
wss_per_k[i] <- kmeans_aux$tot.withinss
}
Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)Warning: Quick-TRANSfer stage steps exceeded maximum (= 3053450)
par(mfrow = c(1, 1))
plot(1:10, wss_per_k,
type = "b",
xlab = "Number of clusters",
ylab = "WSS",
)
Como se puede observar en la gráfica anterior, la suma de cuadrados internos disminuye a medida que aumenta el número de clusters. Sin embargo, a partir de 2 clusters, la disminución de la suma de cuadrados internos es muy pequeña. Por lo tanto, se decide utilizar 2 clusters. Comprobamos que tiene sentido utilizar 2 clusters, ya que conocemos que el dataset es binario.
Generamos el modelo de k-means con 2 clusters.
km_model <- kmeans(df, center = 2, nstar = 20)
Representamos la distribución de los datos en función de los clusters obtenidos.
df$cluster <- factor(km_model$cluster)
plot_ly(df,
x = ~cap.diameter, y = ~stem.height,
z = ~stem.width, color = ~cluster
) %>%
add_markers(size = 1.5)
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Se puede observar que los champiñones de menor tamaño (incluyendo diámetro, altura y anchura) pertenecen al cluster 2 y los de mayor tamaño pertenecen al cluster 1.
A continuación, calcularemos el valor promedio de las variables para cada cluster generado con el modelo de k-means.
grouped_mushroom <- df %>%
group_by(cluster) %>%
summarise(
mean_cap_diameter = mean(cap.diameter),
mean_stem_height = mean(stem.height),
mean_stem_width = mean(stem.width)
)
grouped_mushroom
A partir de este momento, hemos decidido modificar el dataset actual debido a que para aplicar técnicas como “silhouette” o “dendrogram” (para el caso de clustering jerárquico) es necesario que el dataset sea de menor tamaño. Para mantener la proporción de los datos de cada clase, se reducirá el dataset haciendo uso de createDataPartition, manteniendo únicamente un 1% de los datos iniciales para el análisis.
set.seed(42)
split <- createDataPartition(mushroom$class, p = 0.01)
smaller_df <- mushroom[split$Resample1, ]
Comprobamos que la proporción de datos de cada clase se mantiene al hacer la partición.
initial_class_prop <- table(mushroom$class) / nrow(mushroom)
smaller_class_prop <- table(smaller_df$class) / nrow(smaller_df)
print(initial_class_prop)
0 1
0.4450867 0.5549133
print(smaller_class_prop)
0 1
0.4386252 0.5613748
Ya podemos trabajar con el dataset reducido. Comprobamos que el dataset ahora cuenta con 611 ejemplos y 3 atributos.
smaller_df <- smaller_df[, numerical_features]
dim(smaller_df)
[1] 611 3
En primer lugar, visualizaremos la distribución inicial de los datos a través de una gráfica 3D.
plot_ly(smaller_df,
x = ~cap.diameter, y = ~stem.height,
z = ~stem.width
) %>%
add_markers(size = 1.5)
A continuación, vamos a estudiar cuál sería el número óptimo de clusters para el dataset reducido haciendo uso de la medida de bondad interna “silhouette”. Para ello, utilizaremos la función fviz_nbclust de factoextra. Silhouette es una medida que sirve para validar el número de clusters. Se calcula como la diferencia entre la distancia media de un punto a los puntos de su propio cluster y la distancia media de un punto a los puntos de su cluster más cercano.
fviz_nbclust(smaller_df, FUNcluster = kmeans, method = "silhouette")
Según la gráfica, podemos afirmar que el número óptimo de clusters es 2.
Por último, visualizamos la distribución de los datos en función de los clusters obtenidos.
km_sm_model <- kmeans(smaller_df, center = 2, nstart = 20)
cluster <- factor(km_sm_model$cluster)
plot_ly(smaller_df,
x = ~cap.diameter, y = ~stem.height,
z = ~stem.width, color = ~cluster
) %>%
add_markers(size = 1.5)
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Por último, calculamos el valor promedio de las variables para cada cluster generado con el modelo de k-means.
grouped_sm_mushroom <- smaller_df %>%
mutate(cluster = cluster) %>%
group_by(cluster) %>%
summarise(
mean_cap_diameter = mean(cap.diameter),
mean_stem_height = mean(stem.height),
mean_stem_width = mean(stem.width)
)
grouped_sm_mushroom
A continuación, vamos a aplicar el algoritmo de clustering jerárquico a nuestro dataset reducido. Para ello, utilizaremos la función hclust. Primero, calculamos la distancia entre los puntos del dataset.
distance <- dist(smaller_df)
hc_model <- hclust(distance)
Representamos el dendrograma para visualizar la distribución de los datos en función de los clusters obtenidos.
dend_modelo <- as.dendrogram(hc_model)
plot(dend_modelo)
Hasta ahora, hemos obtenido la jerarquía de los datos, pero lo que realmente nos interesa es la clasificación de los datos en función de los clusters. Cortaremos el dendrograma en un punto que nos interese para obtener los clusters. En este caso, hemos decidido cortar el dendrograma en 90 para obtener una visualización del dendograma cortado.
cut <- 0.9
dend_modelo %>%
color_branches(h = cut) %>%
color_labels(h = cut) %>%
plot()
Para obtener el número óptimo de cluster, haremos uso de la medida interna de bondad silhouette. Para ello, utilizaremos la función fviz_nbclust de factoextra.
fviz_nbclust(smaller_df, FUNcluster = hcut, method = "silhouette")
Comprobamos que en este caso, el número óptimo de clusters podría ser 2 o 3, ya que el valor de silhouette es muy similar para ambos casos. En este caso, hemos decidido utilizar 2 clusters para poder comparar posteriormente los resultados con los obtenidos con el algoritmo de k-means.
Calculamos la agrupación del modelo en función del número de clusters que hemos decidido utilizar. Además, calculamos el promedio de los datos de cada cluster para ver si podemos sacar alguna conclusión.
jq_cluster <- cutree(hc_model, k = 2)
grouped_mushroom <- smaller_df %>%
mutate(cluster = jq_cluster) %>%
group_by(cluster) %>%
summarise_all(mean)
grouped_mushroom
Visualizamos la agrupación de los datos en función de los clusters obtenidos.
jq_cluster <- factor(jq_cluster)
plot_ly(smaller_df,
x = ~cap.diameter, y = ~stem.height,
z = ~stem.width,
color = ~jq_cluster
) %>%
add_markers(size = 1.5)
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Warning: minimal value for n is 3, returning requested palette with 3 different levels
Con el objetivo de comparar los resultados obtenidos en los dos algoritmos, vamos a calcular el rendimiento de cada uno de ellos, haciendo uso del accuracy como medida de bondad externa.
En primer lugar, calculamos el accuracy del modelo de k-means. Supondremos que la clase 1 es la clase “e” y la clase 2 es la clase “p”. Para ello, obtenemos las clases reales y las clases predichas, y calculamos el accuracy.
Volvemos a obtener el dataset reducido para poder tener las clases reales.
smaller_df <- mushroom[split$Resample1, ]
real_classes <- ifelse(smaller_df$class == "e", 1, 2)
predicted_classes <- km_sm_model$cluster
predicted_classes <- as.numeric(predicted_classes)
accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)
[1] 0.700491
Hacemos lo mismo con el modelo de clustering jerárquico, pero en este caso, supondremos que la clase 1 es la clase “p” y la clase 2 es la clase “e”.
real_classes <- ifelse(smaller_df$class == "e", 2, 1)
predicted_classes <- as.numeric(jq_cluster)
accuracy <- sum(real_classes == predicted_classes) / length(real_classes)
print(accuracy)
[1] 0.9885434
Tras comparar los resultados obtenidos en los dos algoritmos, podemos afirmar que el modelo de clustering jerárquico ha obtenido un accuracy mayor para este dataset.